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procedure mxinvert (a, n, eps, singular) ; value n,eps; array a; integer n; real eps ; label singular- 
begin integer i, j,k,pivi,pivj,p,ri,ci,rk,cj,iless1 ; real pivot; integer array r,o[1 :nT1 

comment set row and column index vectors ; 

for i:=1 step 1 until n do r[i] : =c[i] : =i ; 

comment find initial pivot; pivi:=pivj:=l ; pivot : =a[1 , 1 ] ; 
f °r i:=1 step 1 until n do for j:=l step 1 until n do 

if abs(a[i, j])>abs(pivot) then begin pivi:=i; pivjr^j; pivot:=a[i, j] end; 
comment start reduction; 
for i:=1 step 1 until n do 

begin ri:=r [pi vi]; r [pivi] : =r[i] ; r[i]:=ri; ci : =c[pivj] ; c[pivj] : =c[i] ; c[i]:=ci; ilessl :=i-1 ; 
if eps > abs(a[ri , ci]) then 

begin print punch(3) ,sameline,digits(3) ,££1?MATRIX SINGULAR?, ££l?i=?,i,££l?PIVOTS FOLLOW? 
for i : =1 step 1 until n do 

print punchC 3) , s amel ine , digi t s ( 3) , ££1?? , r [ i ] , ££s4?? , c [ i ] ; 

goto singular 

end ; 

for j:=1 step 1 until ilessl ,i+1 step 1 until n do 

be e in cj:=c[j]; a[ri,cj]:=a[ri,cj]/pivot end ; 
a[ri,ci]:=l,0/pivot; pivot :=0; 
for k:=1 step 1 until ilessl ,i+1 step 1 until n do 
begin rk:=r[k]; 

for j:=1 step 1 until ilessl ,i+1 step 1 until n do 

begin c j:= c U3r-aC3E-k,cj]:=aCrk,cj3-aCri,cjT*a[rk,ci]; 

if k>i and j>i and abs(a[rk,cj]) >abs(pivot) then 

be g in Pivi:=k; pivj:=j; pivot : =a[rk, cj] end conditional 

end jloop; 

a[rk,ci] :=-aEri,ci]*a[rk,ci] 
end kloop 
end iloop and reduction; 
comment rearrange rows; mxperm(a[ j ,p] ,a[k,p] , j,k,r,c,n,p) ; 

comment rearrange columns; mxperm<a[p , j] , a[p,k] , j,k,c,r ,n,p) 
end mxinvert ; 



LEI 6 



procedure choleski (a, n, fail) ; 

value n; integer n; array a; label fail; 

begin integer i , j,k, ilessl ; real aii; 
for i := 1 step 1 until n do 
begin ilessl := i-1; 

aii:=a[i,i]-sigma(a[i,k]*a[i,k],k,1 ,iless1 ) ; 

if aii < then goto fail; 

aii:=a[i,i]:= sqrt(aii) ; 

for j:= i+1 step 1 until n do 

a[j,i]:= <a[ j,i]-sigma(a[i,k]*a[j,k],k,1 ,iless1 ))/aii 



end 



end choleski 



procedure linv(a,n) ; 

value n; integer n; array a; 

begin integer i,j,k, ilessl ; real aii; 

for i : =1 step 1 until n do 

begin ilessl : a i-1 ; aii : =a[i , i ] ; 

for j:= 1 step 1 until ilessl do 

a[i, j]:= -sigma(a[i,k]*a[k, j],k,j, ilessl )/aii; 
a[i,ij:= 1.0/aii 



end 
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end linv; 



procedure mxmult(a,b,c,m,n,p) ; value m,n,p; integer m,n,p; array a,b,c- 
comment c[m p]:= a[m n] b[n p]; ' 
begin integer i , j , k ; 

for i := 1 step 1 until m do 

for j:= 1 step 1 until p do c[i , j] : =sigma(a[i ,k]*b[k, j] ,k, 1 ,n) 



endmxmult 
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procedure SYMDET(a,n,symdGt,fail) ; value n; integer n; real symdet ; array a; label fail; 
begin integer i, j,k, ilessl ,ci,ii,cj,ij; 

real det ,aii , aki , ai j ; 

det : =1 . ; 

for i:=l step 1 until n do 

begin ilessl :=i-1 ; ci:=c[i]; ii:=i+ci; aii:=a[ii]; 
for k:=1 step 1 until ilessl do 

begin aki :=a[k+ci] ; aii : =aii-aki*aki end; 
if aii<0.0 then goto fail; 
det:=det*aii; 
aii : =a[ii] :=sqrt(aii) ; 
for j:=i+1 step 1 until n do 

begin cj:=c[j] ; ij:=i+cj; aij:=a[ij]; 

for k:=1 step 1 until ilessl do ai j:=ai j-a[k+ci ]*a[k+c j] ; 
a[ij3:=aij/aii 

end j 

end i ; 
symdet : =det 

end SYMDET; 
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procedure SYMSOL(a,b,n) ; value n; integer n; real array a,b; 
begin integer i,j,jlessl ,cj; real bi,bj; 

for j:=1 step 1 until n do 

begin bj:=b[j] ; jlessl:=j-1; cj:.-c[j]; 

for i:=1 step 1 until jlessl do bj :=bj-a[i+cj]*b[i] ; 
b[j] :=bj/a[j+cj3 

end ; 

for i:=n step -1 until 1 do 
begin bi:=b[i]; 

for j:=i+1 step 1 until n do bi : =bi-a[i+c[ j] ]*b[ j] ; 
b[i]:=bi/a[i+c[i]] 

end 

end USYMSOL; 



procedure SYMDET(a,n, symdet , f ail) ; value n; integer n; real syradet ; array a 
begin integer i , j ,ki ,ii,k1 ,kiless1 , ij ,kj ; 

real det , aii , aki , ai j ; 
det:=l .0; ii :=1 ; 
for i : =1 step 1 until n do 

begin aii: =a[ii]; k1:=ii-i+1; kilessl : =ii-1 ; 
" for ki:=k1 step 1 until kilessl do 

begin aki:=a[ki]; aii : =aii-aki*aki end ; 
if aii<0 then goto fail; 

det: =det*aii ; 
aUl=a^ii] :=3qrt(aii) ; 

for j:=iil step 1 until n do 

begin aij:=a[ij3; kj:=ij-i+1 ; 

for ki:=k1 step 1 until kilessl do 

begin ai j : =ai j-a[ki]*a[kj] ; kj:=kj+1 

a[ij]:=aij/aii; ij:=ij+j 

end ; 
ii :=ii+i+l ; 

end ; 
symdet : =det 

end SYMDET; 



procedure SYMS0L( a , b , n) ; value n; integer n; array a,b; 
begin integer i,j,k,ii; real sum; 
ii:=1; 

for i :=1 step 1 until n do 

begin k: =ii-1 ; sum:=b[i]; 

for j:=i-l step -1 until 1 do 

begin sum: =sum-b[ j]*a[k] ; k:=k-l end ; 

b[i] : =sum/a[ii] ; ii:=ii+i+1 

end forward solution; 

k:=ii-1; ii:=ii-n-l; 

for i:=n step -1 until 1 do 

begin sum:=b[i]; 

for j:=i+l step 1 until n do 

begin sum:=sum-b[j]*a[k] ; k:=k+j end ; 

b[i] :=sum/a[ii] ; k:=ii-l; ii:=ii-i 

end back substitution 
end SYMSDL; 



label fail; 



end ; 
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< I 

i 

LL08 ) 

procedure vecjacobi( a, s,n,rho) ; value n,rho; real rho; integer n; array a,s 

comment an adaptation of ACM85 to evaluate the eigenvalues and 

eigenvectors of a real symmetric matrix A[1 :n, 1 :n] . The upper 

triangle of A should be supplied, in vector form, stored by 

columns, in a[1:n (n+1) 2] so that A[i,j] occupies a[i+j 2], 

Alternatively, by symmetry, the lower triangle of A, stored by 

rows may occupy a. At exit the eigenvalues occupy a[1] through 

a[n], with corresponding eigenvectors in the columns of s[1 :n,1 :n] . 

rho is the precision tolerance as used in ACM85, which, 

in practice, should not be less than the relative machine precision. 

On a test matrix of order 20, using a tolerance of 10 -6 vecjacobi 

proved to be 2.5 times faster than ACM85; 

begin integer array c[1:n]; integer i,j,ci,cj,p,q,cp,cq, jlessl , qlessl ,ip,iq; 
switch ss : =main, mainl ; 

real fac,aij,thr,norm1 ,norm2,apq,app,aqq,m,mu, lambda, cost, sint,aip,aiq, 

sip,siq,sincos ,cs45; 

boolean ind; 

cs45:=1 .0/sqrt(2.0) ; 

p:=0; fac:=0.0; 

for i:=l step 1 until n do 

be gi" s[i , i] : =1 .0; c[i]:=p; p:=p+i; 

for j:=i+1 step 1 until n do s[i, j]:=s[j,i] :=0.0; 

end ; 

for j : =2 step 1 until n do 
begin cj:=c[j]; jlessl : =j-1 ; 

for i:=1 step 1 until jlessl do 

begin ai j : =a[i+c j] ; f ac: =f ac+2.0*ai j*ai j end 

end ; 

thr : =norm1 : =sqrt(f ac) ; norm2 : =rho*norm1/n; 
main: thr:=thr/n; 
mainl : ind : =false ; 

for q: =2 step 1 until n do 
begin cq:=c[q]; qlessl :=q-1 ; 

for p:=1 step 1 until qlessl do 
begin apq: =a[p+cq] ; 

if abs(apq) > thr then 
begin cp:=c[p]; ind: = true ; 

app: =a[p+cp] ; aqq: =a[q+cq] ; m:=app-aqq; 
mu : =abs (m) ; 

if mu<rho then cost:=sint:=cs45 

else begin lambda: =sign(m)*apq; mu:=0.5*mu; 

fac:=0. 5/sqrt(lambda*lambda+mu*mu) ; 
cost : =sqrt (0 . 5+mu*f ac) ; 
sint : =lambda*f ac/cost 

end ; 

for i : =1 step 1 until n do 
begin ci : =c[i] ; 

if i < p then begin ip:=d+cp; iq:=i+cq end 
else begin ip:=p+ci; 

iq:=if i>q then q+ci else i+cq 

end ; 

aip:=a[ip]; aiq:=a[iq]; 
sip:=s[i,p] ; siq:=s[i,q]; 
s[i,p] :=cost*sip+sint*siq; 
s[i ,q] : =sint*sip-cost*siq; 
a[ip] : =cost*aip+sint*aiq; 
a[iq] : =sint*aip-cost*aiq 

end i ; 

sincos : =sint*cost ; f ac : =( apq+apq) *sincos ; 

sint : =sint*sint ; cost : =cost*cost ; 

a[p+cp] : =cost*app+sint*aqq+fac; 

a[q+cq] : =sint*app+cost*aqq-f ac ; 

a[p+cq]:=0,0 

end 

end 

end ; 

if ind then goto mainl else if thr>norm2 then goto main; 
for i:=2 step 1 until n do a[i] : =a[i+c[i] ] 
end vecjacobi ; 



real procedure equipol(xbase,y,arg,n,in,li) ; value xbase,arg,m.n,h; real xbase,arg,h; 
begin integer i,j, mlessl; real jh,fi; array f[0:m]; 

if m>n then m:=n; i : rentier ( (arg-xb as e)/h) -m div 2; 
j : - if_ i<G then else if i+m>n then n-m else i ; 
for i:= step 1 until m do f [i] : =y[i+j] ; 
arg: =arg- j*h-xbase ; 
mlessl : =m-l ; 

for i:=0 step 1 until mlessl do 
begin fi:=f [i j ; jh:=h; 

for j:=i+1 step 1 until m do 

begin f[j]:=fi+arg*(f [j]-fi)/jh; jh:=jh+h end ; 
arg : =arg-h 

end ; 
equipol : =f [m] 
end equipol; 



real procedure ait(z,f ,arg,n) ; value arg,n; integer n; real arg; array z,f; 
begin integer i,j,nless1; real fi,zi,u; nlessl : =n-1 ; 
for i:=0 step 1 until nlessl do 

begin fi:=f [i] ; zi:=z[i]; u:=arg-zi; 
for j:=i+1 step 1 until n do 
f[j]:=fi+u*(f[j]-fi>/(z[j3-zi) 

end ; 
ait:=f [n] 

end; 
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?ay y; integer m,n; 
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real procedure equidydx(xbase,y,arg,n,m,h,est) ; value xbase,arg,m,n,h; real xbase,arg,est,h; array y; integer m,n; 
begin integer i,j, mlessl ; real jh,f i ,dif fi ,f jf i ; array f,diff[u:m]; 
if m>n then m:=n; i : =entier< (arg-xbase)/h) -m div L; 
j:= if i<u then else if i+m>n then n-m else T~; 

for i:= step 1 until m do begin f [i] :=y[i+j] ; diff[i]:=0.0 end ; 
ar g : =arg-j *h -xbas e ; 
mlessl : =m-l ; 

for i:=0 step 1 until mlessl do 

begin fi:=f[i]; jh:=h; dif f i : =dif f [i] ; 
for j:=i+l step 1 until m do 
begin f jfi : =f [ j] -f i ; 

diff[j]:=diffi+(f jfi+arg*(diff[ j]-diffi))/jh; 

f[j]:=fi+arg*f jfi/jh; jh:=jh+h 

end ; 
arg:=arg-h 

end ; 
est :=f [m] ; 
equidydx : =dif f [m] 
end equidydx; 



